Proposed realization of itinerant ferromagnetism in optical lattices 
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We propose to realize the itinerant ferromagnetism of two-component cold fermionic atoms in 
the p-orbital bands in optical lattices. The band flatness in the two-dimensional honeycomb lattice 
dramatically amplifies the interaction effect driving the ferromagnetic transition even with a rela- 
tively weak repulsive interaction. This scheme has the advantage that the stability of the system 
can be maintained without suffering decaying to the molecular state as one approaches the Feshbach 
resonance from the side with positive scattering length. Experimental signatures and detections are 
also discussed. 
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Itinerant ferromagnetism(FM) is one of the central 
topics in condensed matter physics [1^7], whose mech- 
anism remains a controversial problem. The major diffi- 
culty is that FM does not have a well-defined weak cou- 
pling picture, i.e., spontaneous spin polarization requires 
interactions strong enough to overcome the kinetic en- 
ergy cost. Furthermore, unlike most superconductors in 
which fermionic excitations are gapped except at nodal 
points which are with zero measure, ferromagnets have 
entire gapless Fermi surfaces. The gapless particle-hole 
excitations dramatically complicates the nature of the 
transition. An important progress in understanding fer- 
romagnetism has been made by Mielke [8] and Tasaki [9] , 
who provide a class of exactly solvable models with FM 
ground states. The key point is the existence of Wannier- 
like orbitals for single particle eigenstates which reduces 
the kinetic energy cost of spin polarization to zero. How- 
ever, these models suffer the stringent conditions of fine- 
tuned long range hopping, which are difficult to realize. 

In spite of its importance, itinerant ferromagnetism has 
not been realized in cold atom systems partly because 
of the requirement of strong repulsion mentioned above. 
It has been proposed [10] to achieve ferromagnetism by 
approaching the Feshbach resonance from the side with 
positive scattering length. However, a problem is that 
before the interaction is strong enough to drive the tran- 
sition, the system is unstable toward to the formation of 
the dimer-molecules [11]. 

On the other hand, orbital physics in optical lattices 
has received a great deal of attention, which gives rises 
a variety of new states of matter with both cold bosons 
and fermions [13-26]. In particular, it has been recently 
shown that the p-orbital band in the honeycomb lat- 
tice, a Pa;,i/-orbital counterpart of graphene, has a natu- 
ral excellent flat band structure. The consequential non- 
perturbative interaction effects (e.g. Wigner crystalliza- 
tion) associated with the band flatness have been inves- 
tigated [17, 18]. Experimentally, the honeycomb lattice 
has been constructed quite some time ago [27]. 

In this article, we propose to realize the ferromagnetic 



state in the p-orbital band of the honeycomb optical lat- 
tice by taking advantage of its band flatness. Due to 
the divergence of the density of states, even weak repul- 
sions can drive the ferromagnetic transition and ensures 
the system stability escaping from the formation of dimer 
molecules, thus overcome this experimental difficulty. We 
show the existence of the exact ferromagnetic ground 
state at the filling level 0.25 < (n) < 0.5, and investi- 
gate the phase boundaries between ferromagnetic, para- 
magnetic, and antiferromagnetic states through the self- 
consistent mean-field method. This research will open 
up a new opportunity to investigate ferromagnetism with 
precise controllability. 

We begin with the free part of the pa;.j,-orbital band 
Hamiltonian in the honeycomb lattice 

r A,i,{T r(7 

where e\^2 = ^^^x + \(iy and 63 = —Cy are the unit 
vectors pointing from an A-site to its three i?-site neigh- 
bors; Pi = {pxfix + PyE-y) ' (« = 1 ^ 3) arc the pro- 
jections of the p-orbitals along the ei direction; rif^.o- is 
the number operator for spin cr; is the chemical poten- 
tial. As illustrated in Ref. [17, 18], its band structure 
contains both flat bands and Dirac cones whose spectra 
are symmetric respect to the zero energy. Two disper- 
sive bands in the middle have two non-equivalent Dirac 
points with a band width of The bottom and top 
bands turn out to be completely flat in the absence of 
the TT-bonding t±, which means that the single particle 
eigenstates of Eq. 1 can be constructed as a set of degen- 
erate localized states. There exists one such eigenstate 
in the bottom band for each hexagon plaquette, whose 
orbital configuration on each site is along the tangential 
direction around the plaquette. The anisotropy of the tr- 
bonding and the destructive interference together forbid 
the particle "leaking" outside, which renders these states 
eigenstates. 

The interaction of the spinful neutral fermions is de- 
scribed within the s-wave scattering approximation, lead- 
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ing to a two-band Hubbard model constructed as 



Hint — U^^(^Tl^^x^'^r,x 



nf^nfy) + Mplx-fPl^iPfyiPfyt + h.c.) (2) 



where U = 



^^a, J Ipxirjl'^d^r (a^ the scattering 
length), J = A = ^, SVa: = |pLa'^a/3P™/3 the 
spin operator at site f in the pj;-orbital, and Spy can be 
defined accordingly. The J7-term is just the usual Hub- 
bard term; the J-term represents Hund's rule physics: 
the on-site repulsion between two fermions in the spin 
triplet states is zero, while that of spin singlet is J; the 
A-term describes the pair scattering process between the 
Px and py-orbitals on the same site. We first consider two 
fermions on the same site to gain some intuition. If they 
are in the spin triplet channel, their orbital wavefunction 
is anti-symmetric, thus the s-wave interaction vanishes. 
If they are in the spin singlet channel, their orbital wave- 
functions are symmetric as Px + Pyi Px ~ Py^ PxPy, 
respectively. The first one has energy U + A — ^ , while 
the later two are degenerate with energy U—A = J=^. 

At low filling factors, (n) < ^, the ground states are 
heavily degenerate. Each fermion occupies one localized 
plaquette state. The touching of two neighboring pla- 
quettes with the same spin costs zero repulsion energy, 
while that of the opposite spin costs energy of Ac- 
cording to the percolation picture of the fiat band fcr- 
romagnetism [8, 9], fermions form disconnected clusters 
composed of touching plaquettes. Let us denote the total 
number of sites N. There are then plaquettes whose 
centers form a triangular lattice. These clusters are la- 
beled from 1 to Tic- The number of plaquettes inside 
the i-th cluster is denoted as m; (i — 1 Uc). Thus 
the particle density is just (n) = X^I^i"^*/-^- Within 
each cluster fermions are fully polarized, but the polar- 
izations of different clusters are uncorrelated. The fer- 
romagnetic order parameter when (n) < ^ is defined as 



M — ^ {Sl^^f)/N, where the square of total spin of the 

system is just S1^g^ — (|)^ X]"=i '^i ' means average 
over all possible configurations of the degenerate ground 
states. As the filling increases beyond the percolation 
threshold, clusters become connected driving the system 
from paramagnetic state to ferromagnetic state. We have 
performed numeric simulations and obtained M v.s. (n) 
as depicted in Fig.l. M becomes nonzero approximately 
beyond the percolation threshold Uc — 0.25, (the small 
deviation is due to finite size effect), which corresponds 
to the site percolation threshold Pc = 0.5 for the triangu- 
lar lattice. As we go beyond the percolation threshold, 
there is an unique infinite cluster with non-zero probabil- 
ity in triangular lattice [28], thus M is just the probabil- 
ity of a site belonging to the infinite cluster and behaves 
as M oc ((n) — Uc)^ with the critical exponent /3 = 
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FIG. 1: The magnetization M v.s. p for < (n) < i simu- 
lated on a triangular lattice with 9000 x 9000 lattice sites(blue 
line) and 300 x 300 lattice sites (red line). Note the strong fi- 
nite size efltect around the percolation threshold (n)c = 1/4. 



[28]. As (n) increases further, M is linear with (n), which 
means that the system is almost fully polarized. 

After the lowest fiat band is fulfilled with spin polar- 
ized fermions, i.e., (n) > i, an exact solution is no longer 
available. We perform a self-consistent mean-field calcu- 
lation. To facilitate the decomposition of the interac- 
tion term of Eq. 2, we can recast it in an equivalent 



form Hint = Ef^Wn + Ef^ 



/^(l)^(l) +^(2) (2), 



where the pseudo-spin operators tP'J'^^ are defined as 



.(2) - 



plt.lPyt,i)'^,i = i(plT,;P^t,; +Pyt,iP2't4)' respectively. 
Here r^'^ measures the preferential occupation of a;-and 
2/-orbitals, / is half of the identity operator in the px 
and Py orbitals. Although this does not make the spin 
SU(2) symmetry explicit, it is convenient for the mean- 
field decoupling if the spin quantization axis is chosen 
along the z-direction. The mean-field Hamiltonian then 
reads, apart from a constant part. 
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The ferromagnetic order parameter is defined as M = 
^-.(/^ — /^) /iV. We enlarge the unit cell to include 
six lattice sites to allow site-dependent magnetic and or- 
bital configurations. We also replace the bare value of 
U with an effective one in the self-consistent calculation 
by taking into account the Kanamori corrections defined 
as Ueff in) = U{p,q;fi) a.t p,q ^ [29, 30] . U{p,q;fi) is 
defined as U{p, q; ~ U/[l — UW{p, q)], and 

H = [E{q) + E{p) - E{k) - Eip + q- k)] , (4) 



where E{k) is the band energy corresponding to Eq.l, 
and the summation over k is taken over unoccupied Bloch 
states. The renormalization procedure ensures that Ueff 
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FIG. 2: A) Filling factor n and magnetization m v.s the 
chemcial potential fi at weak interaction {U — 1). B) n, m 
and antiferrmagnetic moment maf v.s. /i at strong interation 
((7 = 5). 



can be of the same order as the band width at most 
and thus precludes ferromagnetism in the usual Hubbard 
model at low density even for very large bare U [30]. 

As (n) increases beyond i, the flat band is fulfilled 
with polarized particles with /Lt/ty > —1.5. The addi- 
tional particles have to either flip their spin and fill the 
fiat band states costing energy of order U, or they con- 
tinue to fill the second band while remaining spin polar- 
ized. The competition between the kinetic energy and 
the interaction energy determines a critical chemical po- 
tential as a function of U, beyond which the full po- 
larization is no longer favorable. Within the mean field 
theory, we compare the energies of fully polarized ferro- 
magnetic states EferroiU, /i) and that of a paramagnetic 
state Epara{U, fJ-), and then determin /i^. As depicted in 
Fig. 2 A, we show the dependence of (n) and M as a func- 
tion of /z at a weak coupling U /iy = 1. Note that (n) and 
m both jump at ~ —1.15 which is a consequence of 
the fiat band structure. When the minority spins begin to 
populate, the divergent density of states of the fiat band 
suddenly increases its filling, which significantly changes 
the effective potential that the majority spins feel, and 
drive the system to paramagnetic states. In other words, 
if the particle density is fixed between the two values 
before and after the jump, we obtain the phase separa- 
tion state as coexistence of the fully polarized state and 
the paramagnetic state. After the paramagnetic phase is 
reached, the system remains in metallic state at small U 
because of the vanishing of density of states of the Dirac 



cone at half filling. 

The phase diagram becomes richer as U goes stronger 
as depicted in Fig. 2 B with (7/i|| = 5. The fully 
polarized state survives to a larger iic/t\\ ~ 0.1 as ex- 
pected; the Mott-insulating antiferromagnetic state oc- 
curs at half-filling (n) ~ 2 with the Neel order defined 
as Maf — jfJ2ii~y'^l- strong coupling limit at 

half-filling, each site is filled with two particles with spin 
aligned by Hund's rule as a spin-1 complex. The vir- 
tual hopping of particles leads to the antiferromagnetic 
exchange and the long range Neel order in the ground 
state. At intermediate fillings between the fully polarized 
state and antiferromagnetic state, as fi > fic, the quick 
increase of (n) is due to the filling of the minority spins 
in the flat band. Interestingly, an incompressible state 
occurs at (n) = 5/3 with (rt-^) = 1 and (n^) = 2/3, which 
has partial ferromagnetic and partial antiferromagnetic 
order. It also exhibits intricate orbital ordering pattern 
as in the spinless case [17, 18]. We will defer a detailed 
analysis of this state to a future publication. 

Now we discuss the effect of the small 7r-hopping term 
t±. Then the lowest energy band acquires a weak dis- 
persion with band width of order t±. In this case, the 
localized plaquette states are no longer eigenstates of the 
system. However, so long as the interaction strength U 
is much larger than t±, it is not energetically favorable 
to flip the spin and thus the ferromagnetic state is sta- 
ble. As shown in Ref. [18], in realistic systems with the 
sinusoidal optical potential, it is easy to suppress t±/t\\ 
to 10^^ or even 10^'^. We have confirmed the existence 
of ferromagnetic order numerically by using the above 
mean-field theory for t± up to O-lty. 

Another small interaction which was previously left out 
is the p-wave scattering U' . It adds an extra interaction 
term AHint = J2r^'''^fxn-fy, and the value of J in Eq. 
2 changes to J' — — U' . In this case, the flat band 
structure is not destroyed, but an exact solution is no 
long available as long as (n) > 1/6 because different pla- 
quettes touch and cost energy of order of U' even with 
polarized spin. We have also performed self-consistent 
mean field theory to check that the ferromagnetic state 
is still stable at small value of U' . For example, for the 
case of J7/i|| = 5,t/'/i|| = 1, we find that the ferromag- 
netic order survives up to filling factor (n) « 1. 

Now we discuss the finite temperature effect. We only 
consider the case of the filled fiat band whose low en- 
ergy excitations are captured by an effective ferromag- 
netic Heisenberg model defined for plaquette spins [8], 
and leave the study of the general filling to a future pub- 
lication. In our case, the effective model is constructed on 
a triangular lattice as = —J'^Sf Sj with J' ~ j^- 
The mean field calculation based on Eq. 2 shows that the 
critical temperature Tmf is at the same order of zJ'/2 for 
small values of f/iy at {n) = ^ where z = 6 is the co- 
ordination number. T^/ is the temperature scale for the 
onset of the nonzero magnitude of the ferromagnetic or- 
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FIG. 3: The skyrmion configuration of the spin density dis- 
tribution in the confining trap with parameters U/Er = 1, 
t^/ER = 0.24, 9./Er = 0.006, and = = 168. The total 

number of unit colls is A'^ = 24 x 20 (the number of sites is 
2iV). A) {S-,{r)) along the line of {X/a,Y/a = 0), i.e., the 
center row in the system. B) {Sx{r)) pattern along the line 
of {X/a,Y/a = 0). {Sy{r)) shows the similar behavior along 
{X/a = 0,Y/a). 



der parameter, which is roughly one order smaller than 

U. According to the data in Rcf. [33], T„j/ can be esti- 
mated as the order of 1 ~ 10 nK which remains experi- 
mentally accessible. Below T^f, thermal fluctuations of 
the orientation of the ferromagnetic order are described 
by the 0(3) non-linear-cr model, and no true long range 
order exists for T > in 2D. However, the spin corre- 
lation length ^ diverges fast as ^/a w e'^r^f/T jg -j-j^g 
lattice constant) which can easily exceed the small size 
of the optical lattice (e.g. typically L/aK, 50) as T goes 
smaller than T^/ . Thus practically we can still take T^/ 
as the onset temperature scale of ferromagnetism. 

Next we study the eff'ect of the overall harmonic trap- 
ping potential Vir{r) = ^fi^r^ and show that fer- 
romagnetism remains robust through a self-consistent 
Bogoliubov-de Gennes calculation in real space. In ex- 
periments, the particle number of each component is sep- 
arately conserved. Furthermore, when the system has 
SU{2) symmetry as in Eq. 2, the spin configuration 
prefers the skyrmion configuration rather than the Ising 
domain wall [36] to lower the energy. The skyrmion- 
type real space magnetization distributions of {Sz{r)) 
and {Sx{r)) are depicted in Fig. 3 A and B, respectively, 
with the parameters U/t\\ = 1; fi = 0.006£'ii {Er the 
recoil energy); t\\/ER = 0.24 as calculated in Ref. [19]. 

The characteristic trap length scale is I = ^ ^ « 4.1a 

where a is the lattice constant. Wc take particle num- 
bers of spin t and J, components as = = 168 which 
can be achieved by introducing two different chemical 
potentials for /i^ = —1.3897 and /i^ = —1.3903. The 
spin ], particles concentrate in the center region, whereas 
most spin f particles live outside and form a ring. The 
spin distribution smoothly varies from 4- to f from the 
central region to the outside due to the spin SU{2) sym- 
metry. Experimentally, the shell structure can be re- 
solved in situ by tomographic radio- frequency spectro- 
scopic method [34] followed by a phase-contrast imaging 



technique [35] which can distinguish the shell structures 
with two different components. 

In summary, we have proposed a scheme to realize the 
itinerant ferromagnetism by using the flat band struc- 
ture in the p-orbital system in the honeycomb lattice. 
The band flatness facilitates the stabilization of ferro- 
magnetism even at weak repulsive interactions. Anti- 
ferromagnetic Neel ordered states are found at large re- 
pulsive interactions at half-filling. This work potentially 
provides a way to realize spin transport and even spin- 
tronics in the future research with cold atoms. 

S. Z. is supported by NSF-DMR-03-50842. C. W. 
thanks J. Hirsch and L. Sham for helpful discussions. C. 
W. and H. H. Hung are supported by the Sloan Research 
Foundation, ARO-W911NF0810291 and NSF-DMR-03- 
42832. 

Note added After this work was completed, there ap- 
peared independent works by Wang et al. [37] and Bcrd- 
nikov et al. [36] on the ferromagnetism on cold atoms. 
Also evidence of ferromagnetism has been observed in a 
recent experiment [38], which, however, is currently un- 
der intensive debating . The experiment result could be 
understood as a strongly correlated paramagnetic state 
with intensive ferromagnetic fluctuations [39]. 
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